Astronomy & Astrophysics manuscript no. secular 


©ESO 2013 


Januai-y 20, 2013 





Secular dynamics of planetesimals in tight binary systems: 

Application to 7-Cephei 

C.A. Giuppone^ A.M. Leiva^'^, J. Correa-Otto'-^, and C. Beauge^'^ 

' Observatorio Astronomico, Universidad Nacional de Cordoba, Laprida 854, (X5000BGR) Cordoba, Argentina 

e-mail: cristian@oac.uncor.edu 
- Instituto de Astronomfa Teorica y Experimental, Laprida 854, (X5000BGR) Cordoba, Argentina 

Preprint online version: January 20, 2013 

^ ABSTRACT 

. Context. The secular dynamics of small planetesimals in tight binary systems play a fundamental role in establishing the possibility of 

accretional collisions in such extreme cases. The most important secular parameters are the forced eccentricity and secular frequency, 
' which depend on the initial conditions of the particles, as well as on the mass and orbital parameters of the secondary star. 

Aims. We construct a second-order theory (with respect to the masses) for the planar secular motion of small planetasimals and deduce 
new expressions for the forced eccentricity and secular frequency. We also reanalyze the radial velocity data available for y-Cephei 
' and present a series of orbital solutions leading to residuals compatible with the best fits. Finally, we discuss how different orbital 

configurations for y-Cephei may affect the dynamics of small bodies in circunmstellar motion. 

Methods. The secular theory is constructed using a Lie series perturbation scheme restricted to second order in the small parameter. 
■ The orbital fits were analyzed is done with a minimization code that employs a genetic algorithm for a preliminary solution plus a 

' simulated annealing for the fine tuning. 

• ' Results. For y-Cephei, we find that the classical first-order expressions for the secular frequency and forced eccentricity lead to large 

i-C ' inaccuracies ~ 50% for semimajor axes larger than one tenth the orbital separation between the stellar components. Low eccentricities 

and/or masses reduce the importance of the second-order terms. The dynamics of small planetesimals only show a weak dependence 
with the orbital fits of the stellar components, and the same result is found including the effects of a nonlinear gas drag. Thus, the 
^ ' possibility of planetary formation in this binary system largely appears insensitive to the orbital fits adopted for the stellar components, 

^ , and any future alterations in the system parameters (due to new observations) should not change this picture. Finally, we show that 

^ . planetesimals migrating because of gas drag may be trapped in mean-motion resonances with the binary, even though the migration 

is divergent. 

Key words, planets and satellites: formation - stars: binaries: close: y-Cephei - methods: data analysis - methods: analytical - planets 
^ and satellites: dynamical evolution and stability. 

m 

1. Introduction The y-Cephei system then constitutes a paradigm. It may be 

argued that if we are able to comprehend planetary formation in 

Although it is believed that approximately half of all stars belong such an extreme environment, we would have taken large steps 

O to multiple stellar systems (e.g. Duquennoy and Mayor 1991), towards a global understanding of planetary formation in any 

^-H ~ 90% of exoplanets are associated with single stars (Zsom et al. other system. It is then no surprise that this planetary system has 

, 2010). It is not yet clear whether this discrepancy is solely due caught the attention of several researchers over the past decade, 

^ ■ to observational bias, or if the process of planetai-y formation and many dynamical and collisional studies may be found in 

• ' may be seriously impaired even in veiy wide binary systems, the literature (e.g. Thebault et al. 2004, 2006, Haghighipour 

^ Curiously, however, a few exoplanets have also been detected in 2006, Veiiier and Evans 2006, Tsukamoto and Makino 2007, 

H very tight binary stellar systems, where the gravitational pertur- Paardekooper et al. 2008, Xie and Zhou 2008, Kley and Nelson 

- . . bations of the secondary component are so large that accretional 2008, Paardekooper and Leinhardt 2010). 
collisions among small planetesimals are extremely difficult. 

Perhaps the most extreme case is y-Cephei, a binary stellar In Beauge et al. (2010) we showed that the dynamics of the 
system whose most recent orbital determination (Neuhauser et gas disk plays an important role in determining the evolution of 
al. 2007) shows a secondary component of mass rriB ~ 0.4Mq small planetesimals. We found that, although an apsidal preces- 
orbiting a principal star niA ~ 1 .4Mq in an ellipse with semi- sion of the gas elements may play a disruptive role, especially 
major axis ob ~ 20.2 AU and eccentricity eg ~ 0.41. Although in the inner parts of the disk, the combined effects of a non-zero 
both stars have a minimum mutual distance of only ~ 12 AU, a precession rate plus a high forced eccentricity in the disk may 
giant planet has been detected at ~ 2 AU from niA (Hatzes et al. in fact lower the relative velocity of solid bodies in the outer 
2003). So far, all attempts to understand the accretional history regions. We thus envisioned a possible scenario in which the 
of this extrasolar planet have been unsuccessful, and it is diffi- growth from kilometer-size planetesimals to ~ 100 km plane- 
cult to visualize an scenario under which such a massive Jovian tary embryos could initially take place near the truncation radius 
planet could form through accretional collisions from a primor- of the gas disk. As the embryos spiral down towards a lower 
dial planetesimal swarm. semimajor axis, subsequent collisions could then lead to larger 
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embryos, finally forming a giant planet core hopefully near the 
present location of the currently detected Jovian planet. 

To test this idea, we performed preliminary N-body simula- 
tions of the dynamical and collisional evolution of planetesimal 
swarms in the outer regions of the gas disk. However, it was soon 
realized that the dynamics of the solid bodies in this region is ex- 
tremely complex. First, we found that mean-motion resonances 
with the secondary star have significant effects, even though they 
are of very high order (e.g. 12/1, 11/1, 10/1). Second, the prox- 
imity to the secondary star also affects the secular dynamics, 
and the classical analytical models (e.g. Heppenheimer 1978) 
widely used in these problems fail to reproduce the correct or- 
bital variations. The imprecision does not lie in the reduced ex- 
pression adopted for disturbing potential, but in the averaging 
process used to eliminate short-period terms. Similar to the ir- 
regular satellites around the outer planets of our solar system 
(e.g. Cuk and Bums 2004, Beauge et al. 2006), higher order sec- 
ular effects from the interaction of short-period terms (includ- 
ing evection) must be considered to represent the dynamics of 
planetesimals in close binary systems. However, contrary to the 
satellite problem, here the perturber lies in a high-eccentricity 
orbit, and classical high-order models cannot be applied directly 
(Correa Otto et al. 2010). 

The purpose of this work is to present a general description 
of the secular dynamics of small planetesimals in circumstellar 
motion around the more massive star of a generic binary stellar 
system. We assume that all bodies share the same orbital plane. 
Although it is known that even moderate mutual inclinations can 
have significant effects in the accretional evolution of a planetes- 
imal swarm (e.g. Marzari et al 2009a, Xie and Zhou 2009, Xie et 
al 2010, Fragner et al 201 1), here we concentrate on the planar 
case and leave the extension to 3D for a future work. Finally, al- 
though our model will be generic, we apply the results to the par- 
ticular case of y-Cephei where we analyze how the uncertainties 
in the orbital fits of the secondary stellar companion may affect 
the evolution of a planetesimal swarm. 

The paper is organized as follows. Section 2 presents our 
second-order perturbation model and analytical approximations 
for both the forced eccentricity and secular frequency. Since we 
focus our attention on y-Cephei, Section 3 discusses the orbital 
parameters determined for the two stellar components and their 
precision. The secular dynamics of individual planetesimals, un- 
der the additional effects of a nonlinear gas drag, is analyzed 
in Section 4. We also present an example of resonance trapping 
obtained with divergent migration. Finally, discussions close the 
paper in Section 5. 



2. Secular dynamics 

Let us assume a small planetesimal of mass m in circumstellar 
motion around a star of mass niA, which is in turn part of a binary 
system with a smaller component of mass mb- Let be the m^- 
centric semimajor axis of niB and cb its orbital eccentricity. We 
further assume that all motion occurs in a plane. 

Neglecting the gravitational effects of m on both stellar bod- 
ies, the orbit of nis will be a fixed ellipse, while the motion of the 
small planetesimal will be perturbed by the gravitational effects 
stemming from the secondary component. Thus, in our dynam- 
ical system, nis will play the role of the perturber, while m will 
be the perturbed mass. 



2.1. The first-order secular model 

Outside any significant mean-motion resonance, the orbital evo- 
lution of m will be dominated by the secular perturbations, and 
the short-period terms (associated to the mean longitudes) can 
be eliminated by a perturbation technique known as averag- 
ing. The expression for the secular disturbing function R usu- 
ally employed for these studies was originally developed by 
Heppenheimer (1978) which, except for constant terms, is given 
by 



QmB 



8(l-4)3/2fl 



aecB 



2flB(l -e\) 



cos {m — tub) 



(1) 



(see Terquem and Papaloizou 2002), where Q is the gravitational 
constant, a is the m/i-centric semimajor axis of the planetesimal, 
e its eccentricity, and cr its longitude of pericenter. The angle tub 
denotes the longitude of pericenter of the orbit of m^, assumed 
constant. 

Expression ([T]i is constructed from Kaula's (1962) expansion 
of the disturbing potential, truncated to second-order expansion 
in the eccentricity of the perturbed body, and performing a first- 
order "scissors" averaging (with respect to the masses) in the 
mean longitudes. We refer to the resulting expressions as ?l first- 
order model for the secular dynamics. 

Since R does not depend explicitly on the mean longitude 
A of the planetesimal, its semimajor axis is constant and equal 
to the proper value a* . Consequently, the secular system is re- 
duced to a single degree of freedom, and the differential equa- 
tions governing the regular variables k - e cos {m - tub) and 
/i = e sin {m - tub) can be written as 



dk _ 
dt 

where 



-gh 



dh 
dt 



■ 8(k - ef), 
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(2) 

(3) 
(4) 



Given arbitrary initial conditions (ko, ho), these equations ad- 
mit periodic solutions of the form 



k(t) = ej,cos(gt + (po)-\-ef 
h{t) = epSinigt + (f>o). 



(5) 



where g is the secular frequency, e^, = {ko - e f)^ + is usually 
known as the proper (or free) eccentricity, and the phase angle 
is given by the expression tan(^o - h()/{ko - ef). The constant 
term e/ is known as the forced eccentricity and is only present 
in systems with an eccentric perturber Adopting fixed values for 
aB and es, equation (|4]l implies that e/ is a linear function of 
the proper semimajor axis (e/ ~ a*) while the secular frequency 
scales as ^ ~ a*^^^. 



2.2. Numerical simulations 

Our first task is to assess the accuracy of the secular solutions 
(|5]l corresponding to the first-order model ([!]• Two quantities 
we particularly wish to test are e j and g. The forced eccentricity 
is crucial in determining the equilibrium eccentricity of plan- 
etesimals under the effects of gas drag from the protoplanetary 
nebula. Although any secular oscillatory motion is expected to 
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Fig. 1. Forced eccentricity (top) and secular frequency (bottom), 
as function of the proper semimajor axis, calculated by three dif- 
ferent methods: filtered exact numerical simulations (filled black 
circles), semi-analytical first-order averaging of the exact dis- 
turbing function (dashed lines), and the classical analytical first- 
order secular model using Eqs. |3]and|4](continuous lines). 



be damped in a gas-rich scenario, the magnitude of g is impor- 
tant for establishing the validity of the averaging process of the 
disturbing function. 

For our computations, we assume a generic binary system 
with mass ratio between the components of mslniA - 0.4 and 
eccentricity cb - 0.36. This value is similar to the best-fit solu- 
tion found by Hatzes et al. (2003) for y-Cephei. 

Figure [T] shows three different calculations of the forced ec- 
centricity (top frame) and the secular frequency (bottom frame). 
The value of ey appears to grow linearly with the proper semi- 
major axis, reaching values of ~ 0.1 for a* ~ Q.Iqb- The numer- 
ical calculations were obtained from a long-term integration of 
the exact equations of motion, after an online application of a 
low-pass FIR (finite impulse response) filter (e.g. Carpino et al. 
1987). 

A digital filter is a numerical tool that eliminates certain 
frequencies from an input signal. For example, given a certain 
time series (e.g. eccentricity as function of time) and a pass fre- 
quency Vpass, applying a low-pass filter signal will yield an out- 
put that maintains all the periodic variations with frequencies 
y < Vpass while eliminating the rest. Digital filters are a common 
tool for constructing synthetic theories of long-term asteroid dy- 
namics (e.g. Knezeviic and Milani 2000) and planetary dynam- 
ics (e.g. Michtchenko and Ferraz-Mello 2001), and constitute 
a useful alternative to analytical perturbation theories when the 
Hamiltonian function is very complex. 

For the present work, the parameters of the digital filter 
were chosen to eliminate all periodic variations with up to eight 
times the orbital period of the binary component. In a dynami- 



cal system displaying regular motion, the application of the fil- 
ter is equivalent to a full (i.e. infinite-order) averaging of the 
Hamiltonian. The region located beyond a* ~ Q.Iob shows dy- 
namical instabilities that complicate the determination of the 
secular solution. 

A comparison between the numerical and the analytical val- 
ues shows significant differences. Although both methods yield 
similar results for low values of the semimajor axis, the exact 
secular frequencies are systematically underestimated by the an- 
alytical model, leading to differences of almost a factor of two 
for a* /qb ~ 0.24. This limitation in the classical estimation of 
g was first noticed by Thebault et al. (2006), who presented an 
empirical functional correction term to expression (|3). This cor- 
rection reduced the discrepancy to values of around 5% in the 
same range of a*. 

Perhaps more important is that the value of the forced eccen- 
tricity also shows significant differences. While the analytical 
model predicts a mono tonic increase in ey as function of a*, the 
real value appears to reach a plateau around a* /ah - 0.17 (corre- 
sponding to ey ^ 0.07) and decrease for larger radial distances. 
The scatter in the numerical values of both ey and g in this outer 
region stems from the action of high-order mean-motion reso- 
nances between the massless body and the binary star niB- 

At first hand, it seems natural to believe that the limitations 
of expressions ©-(Ell are due to the truncation of the disturbing 
function to second-order terms in the eccentricities and/or third- 
order terms in the ratio a/oB- However, this is not the case. In 
Figure [T] we have also plotted the same quantities determined 
using a semi-analytical model for the disturbing function. This 
expression is calculated directly as 



■2k 



1 



tbI 



■ cos 4>jl'^dAB, 



(6) 



where r and Tb are the position vectors of m and niB, respectively, 
r and are their absolute values, A and Ab are the mean longi- 
tudes, and (p is the instantaneous angular distance between both 
bodies. The integrand is the exact expression for the disturbing 
function with no approximations, and the double integral is per- 
formed numerically. From this expression the value of ey can be 
estimated from the minimum value attained by {R) in the line 
segment (m - tuq) - 0, while the secular frequency is given by 

^ d{G-L) ^ ' 

at the same point. Here (G-L) ^ V^'Waa (e^/2) is the modified 
Delaunay canonical momenta conjugate to the longitude of the 
pericenter Expression (|7) may also be evaluated numerically for 
any initial condition. 

This type of semi-analytical model has proved to be a pow- 
erful tool for mapping the phase space of complex dynami- 
cal systems, especially in the high-eccentricity regime where 
analytical approximations for the Hamiltonian are not avail- 
able (e.g. Michtchenko and Malhotra 2004, Michtchenko et al. 
2006). Formally, it is equivalent to a first-order averaging (in the 
masses) of the exact Hamiltonian function. 

Figure [T] shows the values of both ej and g determined with 
this semi-analytical approach. Although the value of the secu- 
lar frequency shows a significant improvement over the analyt- 
ical estimation, there is still a discrepancy with the exact value. 
This is even more noticeable in the forced eccentricity, where 
there is practically no difference with the value determined from 
equation (O. Consequently, it appears that the limitations of the 
analytical model are not due primarily to the truncation of the 
disturbing function. 
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2.3. A second-order secular model 

Since the errors in the estimation of both the forced eccentric- 
ity ef and the secular frequency g do not come from the Hmi- 
tations in the adopted disturbing function, their origin must lie 
in the construction of the averaged solution itself. As mentioned 
previously, Heppenheimer's (1978) expressions are a first-order 
model with respect to the perturbing mass. Here, we extend the 
calculations to the second order 

One of the most widely used perturbation techniques is the 
so-called Hori's averaging process (Hori 1966, see also Ferraz- 
Mello 2007), which employs Lie-type canonical transformations 
to eliminate the dependence of the Hamiltonian with respect to 
a given set of variables. The new Hamiltonian function is given 
by a power series in the small parameter (e.g. perturbing mass). 

Since we adopt a Hamiltonian formulation, we first need 
to introduce canonical variables. We have chosen the modified 
Delaunay variables (L, A,G - L, A, As, tu), where the canonical 
momenta are given in terms of the orbital elements, by 



G - L = L( Vl - e2 _ 1) 



(8) 



and A is the canonical conjugate of the mean longitude of the 
perturbing mass (i.e. Ab). This third degree of freedom appears 
when passing to the extended phase space to eliminate the non- 
autonomous character of the perturbation. 

The full Hamiltonian function governing the dynamics of the 
planetesimal m is given by 



F(L, A,G-L, A, Ab, ct) 



+ n^A - R 



(9) 



where hb is the mean-motion of the perturbing mass niB and R 
the disturbing function. We can express this Hamiltonian in a 
form adequate for perturbation theory: F - Fq + eF\, where 



Fo = 



aB 
"ir-rel 



(10) 



ruB 



and E - QniBluB is a small parameter that serves as a guide of 
the relative magnitudes between the perturbation term F\ and 
the unperturbed integrable Hamiltonian Fq. 

For the disturbing function, we adopt a Legendre expansion, 
truncated to fourth order in the ratio qIub', in other words, we 
approximate the perturbation by 



f,(cos ^ 



(11) 



where Pi(cos(p) is the Legendre polynomial of degree 
Switching from a power series in cos to a harmonic decom- 
position in (p and transforming them to orbital elements, we can 
obtain a truncated expansion of the disturbing function leading 
to 



Fi ^ ^ ^We'e^g cos {kM + 1Mb - sm) 



(12) 



iJ.S=0kJ=-O3 



where M and Mb are the mean longitudes of both bodies, and 
Dij,ic,i be obtained in terms of the Hansen coefficients (see 
Beauge and Michtchenko 2003). 

Having an explicit expression for F[ in mean variables, we 
may now apply Hori's method. The idea is to search for a Lie- 
type canonical transformation B - eB\ +b^B2+. . . to a new set of 



variables (L*, A*, (G - L)* , A*, A*^, m*) such that the transformed 
Hamiltonian F* is independent of A* and A*^. Up to second order 
in the small parameter, the new Hamiltonian function may be 
written as 



F*((G - L)\ Am*; L", A") ^Fl+ eF\ + s^F^ 



(13) 



where Act* - vj* - tub. The different orders in expression (ITJt 
are given by 



Fq{L\A*) 
{Fi). 



A,Ab 



(14) 



F2 = :^{{{Fi+F\),B^]),,,, 



where {) is the Poisson bracket, {)a,Ab denotes the averaging 
with respect to both mean longitudes (keeping all other variables 
fixed), and B\ is the first-order generating function of Hori's 
method. In terms of the adopted expansion for the disturbing 
function (fT2l i. it is given by 



Bi 



i,j,s,k,l 



F>i,m 

kn + Ihb 



e*'e*J sin (kM* + lM*g - sm*). 



(15) 



where the function must be evaluated in the new variables. 

The construction of the new secular Hamiltonian F*((G - 
L)* , Am* ; L* , A*) is cumbersome, although fairly straightfor- 
ward when using an algebraic manipulator Fortunately, it will 
not be necessary to write an explicit expression here. Let it suf- 
fice to say that F* constitutes a second-order model of the secu- 
lar system and a single degree of freedom system in variables 
((G - L)*,AziT*). Employing the inverse transformation from 
Delaunay variables to orbital elements, we can also obtain an 
expression for F*{e*, Am*; a*) in terms of the mean eccentricity 
e* and the proper semimajor axis a*. Since the latter orbital el- 
ement is constant, it appears in the Hamiltonian as an external 
parameter 

Finally, after solving the secular system and obtaining both 
e* and Act* as functions of time, we may invoke the inverse Hori 
transformation to obtain the short-period variations of the origi- 
nal osculating variables. For the eccentricity, this yields 



e\t) 



.2 2e dBi 
e (f) H T, — 

L* dm* 



(16) 



Because Bi explicitly depends on the mean longitudes, the sec- 
ond term models the short-period variations in the eccentricity, 
while the first term (e*^(f)) gives the main secular contributions. 
Since the eccentricity is a positively defined function, the mag- 
nitude of the second term also specifies the minimum mean ec- 
centricity e* of the secular system for any given proper semima- 
jor axis a*. At the same time, it also gives the averaged semi- 
amplitude of the short-period variations Ae in the same orbital 
element. 

Figure |2] shows an application of our second-order model to 
the same generic binary system as was discussed in Figure [1] 
The plot shows the forced eccentricity, as a function of the ratio 
a* Iub, calculated with three different methods. Recall that this 
model predicts a linear increase of e/ with the semimajor axis. 
Finally, the value of the forced eccentricity determined with our 
second-order Hamiltonian F* is shown as a continuous curve. 
The agreement with the numerical data is very good, and the 
saturation in the value of e/ is reproduced quite well. Since we 
have avoided all small denominators in the generating function 
Bi, the model curve is smooth and shows no indication of the 
effects of mean-motion resonances. 
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Fig. 2. Forced eccentricity as function of the proper semimajor 
axis, calculated by three different methods: filtered exact numer- 
ical simulations (filled black circles), first-order analytical model 
(dashed lines), and the new second-order secular model (contin- 
uous lines). 



2.4. Extending the Thebault et al. (2006) approximation 

As mentioned before, although the second-order model leads to 
significant improvement in the secular solution, as well as al- 
lowing the magnitude of the short-period orbital variations to 
be modeled, it is much too complex to constitute a workable 
model. For this reason, we wondered whether the empiric cor- 
rection term introduced by Thebault et al. (2006) for the secular 
frequency could be extended to reproduce both the forced ec- 
centricity and the short-period variations. Of course it is not ex- 
pected to yield the exact same results, but if the errors are not 
significant, such an empirical second-order approximation could 
constitute a simple quantitative analytical model. 

Following the same approach as Thebault et al. (2006), we 
use e /q and go to denote the first-order expressions for the forced 
eccentricity and secular frequency, and reserve e f and g for the 
second-order values. The idea then is to write e/ = ^foi^ +s6ef) 
(and a similar equation for g), and attempt to model the correc- 
tion terms 6e f and 6g. After several tests and multivariate linear 
regressions, we find that the expressions 



^0 
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(17) 
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Fig. 3. Variation in the forced eccentricity (top) and secular fre- 
quency (bottom), in terms of the proper semimajor axis, for three 
values of the binary eccentricity eg. As before, filled circles 
present results from filtered exact numerical simulations, dashed 
lines correspond to the first-order analytical model, while the 
empirical solutions ( fTTl l are shown in continuous lines. 



Finally, the semi-amplitude of the short-period variations in 
eccentricity can also be empirically modeled according to the 
expression 



\mAj\aBj 



(1 



»2 \6- 



(19) 



Figure[3]once again compares the estimated values of e/ and 
g, this time for three different values of the eccentricity eg of the 
binary component. Given the simplicity of these equations, the 
agreement with the numerical results is surprisingly good. 



agree with the complete second-order model very closely. There 
are some slight differences in g with respect to the original for- 
mula introduced by Thebault et al. (2006) but they are minor and 
not very significant. Finally, the expressions for e jq and go are 
those given in Q and (|4|. 

In terms of ( fTTl l the secular Hamiltonian may be approxi- 
mated well by 



n a g 



(18) 



where k* = e* cos (Am*) and h* - e* sin (Am*) are the new reg- 
ular secular variables. 



3. The y-Cephei binary system 

After specifying the basic ingredients of our second-order dy- 
namical model, we attempt to apply it to y-Cephei. As men- 
tioned in the introduction, this is probably the best-studied tight 
binary system with a known planetary body. Since the main sec- 
ular parameters g and e / strongly depend on the stellar masses 
and orbital elements of the secondary star, we begin our discus- 
sion by reviewing the accuracy of these parameters. 

Throughout this work we refer to the more massive stellar 
component by y-Cephei-A, while y-Cephei-B is used to identify 
the less massive star The giant planet orbiting y-Cephei-A is 
called y-Cephei-b. The masses of each body are denoted by m^, 
niB, and irip, in that order 
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3.1. History and radial velocity data 

Several years before the discovery of the first planetary body 
around a main sequence star (Mayor and Queloz 1995), 
Campbell et al. (1988) suggested the presence of a Jupiter-mass 
object in a 2.7 yr orbit around y-Cephei-A. The authors, how- 
ever, remained cautious about claiming a true planetary detec- 
tion, since the observed periodic variations in radial velocity 
(RV) were at the very limit of the instrumental resolution. To 
complicate the problem even further, the variations in RV at- 
tributed to the Jovian planet, with a semi-amplitude of only about 
25 m/s, were superimposed on a much larger variation caused by 
a previously unnoticed stellar companion with a much longer or- 
bital period. 

The planetary interpretation was questioned later when 
changes in the chromospheric activity were observed with sim- 
ilar period (Walker et al. 1992). Thus, it was proposed that the 
observed changes in RV were spurious and probably only due to 
changes in the spectral line profiles caused by surface inhomo- 
geneities (spots). 

The existence of a binary component (i.e. y-Cephei-B) was 
only reevaluated several years later, when Griffin et al. (2002) 
combined several historical sources of radial velocities that in- 
clude epochs from 1896 to 1980. This data set consisted of 88 
RV observations, although many of them did not contain proper 
uncertainties, and a gap of some 50 years was present in the data 
set. Even so, the authors proposed a secondary stellar mass in 
the system with an orbital period of f ~ 66 yrs. 

The presence of a third body, this time a planet around y- 
Cephei-A, was only confirmed by Hatzes et al. (2003) after in- 
corporating new high-precision velocity observations from the 
McDonald Observatory. They show convincingly that the 2.5 
yr variation in radial velocity was coherent in phase and ampli- 
tude throughout the entire 20 yr interval, as would be expected 
for Keplerian motion, and that no changes were observed in the 
spectral-line bisectors. 

More recently, Torres (2007) has again analyzed the histor- 
ical sources of radial velocities using the extensive Harvard- 
Smithsonian Center for Astrophysics (CfA) database consisting 
of ~ 250, 000 spectra. Torres pointed out that some of the histor- 
ical radial velocities showed large internal discrepancies when 
compared with other data taken at similar times and were con- 
sequently not reliable. The author constructed a reliable data set 
consisting of 30 RV observations. The complete sets of radial 
velocities (four sets by Hatzes et al. (2003) and one by Torres 
(2007)) are shown in Figure |4] where the errors bars indicate the 
uncertainty on each numerical value. The difference in precision 
is remarkable, showing how the incorporation of modem tech- 
niques in RV measurements lead to the detection of the planetary 
mass. This increase in precision also allowed the mean anomaly 
M and longitude of pericenter m of the binary component to be 
accurately defined. 

Independently and without attempting to identify any plan- 
etary body, Gontcharov et al. (2000) studied the mass and or- 
bital parameters of the binary system using astrometric observa- 
tions from several sources. They obtained an orbital period of 
~ 45 yr and a total mass of 3Mq. Unfortunately, the individual 
masses were not specified. In his work, Torres (2007) also com- 
bined a total of 140 astrometric measurements obtained between 
the years 1898 and 1995. Because of the relatively short time 
span of these observations compared to the binary orbital pe- 
riod, Torres noted that there is a high probability that part of the 
orbital motion of the binary has been absorbed into the proper 
motion components reported by Hipparcos. The astrometric in- 




4000 
_ 2000 

CO 







-20000 


-10000 
JD [-2,440,000] 




10000 






1 


1 


1 




' 1 

m = 1.59 M- 

A 






1 


1 


1 




1 






-20000 


-10000 







10000 



JD [-2,440,000] 

Fig. 4. The five sets of RV data used for the orbital fit of y- 
Cephei-B. The four datasets from Hatzes et al. (2003) are shown 
in black, while the dataset from CfA by Torres (2007) is shown 
in gray. The error bars correspond to the observational uncertain- 
ties given by the authors. Two orbital solutions are shown, one 
corresponding to a larger orbit for the binary (top) while the bot- 
tom panel represents a more compact configuration. Each plot 
also assumes a different value for m^. 

Table 1. Published masses and orbital parameters for y-Cephei 





niA 


niB 


ag 


P 


eg 


■njg 






[Mo] 


[AU] 


[years] 




[deg] 


(1) 


niji + nig = 3 




-45 






(2) 


1.59 


0.34 


18.50 


56.81 


0.36 


158.8 


(3) 


1.18 


0.36 


19.02 


66.8 


0.41 


160.9 


(4) 


1.40 


041 


20.18 


67.5 


0.41 


161.0 



(1) Gontcharov et al. 2000, (2) Hatzes et al. 2003, (3) Torres 2007, (4) 
Nehusauer et al. 2007. 



formation from both these papers is quite different (compare Fig. 
6 in Torres 2007 with Fig. 5 in Gontcharov et al. 2000) and the 
discrepancy in the data prior to 1979 prevent us from using any 
astrometric data in our analysis of the orbital determination. 

Although the high-precision RV measurements from Hatzes 
et al. (2003) give very good definitions of the mass and orbital 
parameters of the planetary body, the orbit of binary itself is far 
from established. In Table [1] we summarize the results of four 
different best fits: not only are there noticeable differences in the 
semimajor axis, but the stellar masses also show large discrep- 
ancies. 

The brightness of y-Cephei has made it an easy target for 
spectroscopic studies to determine the effective temperature of 
the star. This parameter, together with the absolute magnitude. 
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are used to determine its mass. However, the effective tempera- 
ture for the star varies from 4300 to 5 100 K (Torres 2007) yield- 
ing a wide variety of possible solutions. 

3.2. Possible orbital solutions for y-Cephei 

Together with the RV data, Figure |4] also shows two different 
orbital fits for the binary system, each constructed for differ- 
ent niA. The top frame adopts the value given by Torres (2007), 
which is significantly lower than the one employed by Hatzes et 
al. (2003), shown here in the bottom graph. We recall that the 
higher mass is usually used for dynamical studies in this system. 

From the raw RV data, we redetermined the best fits for each 
value of niA. To this end we used the PISA code (Pikaia ge- 
netic algorithm + simulated annealing) that we developed for 
our studies of resonant exoplanetary systems (e.g. Giuppone et 
al. 2009). The minimization procedure implies a determination 
of ten free parameters: five for orbital parameters and five for 
the RV offsets. We neglected the presence of the planetary body, 
since it does not introduce any significant effect on the orbital 
calculation for the binary system. The values of ms, as, and es 
obtained for each solution are 

niA = I.I8M0 : niB = O.3IM0 ae = 18.84AU eg = 0.41(20) 
niA = 1.59Mo : iub = O.37M0 ub = 20.70AU eg = 0.41. 

As mentioned previously, the mean anomaly and the longitude of 
the pericenter are well specified, and show little change in both 
fits. We assumed an edge-on configuration. This is an important 
restriction, since any inclination of the orbital plane of the bi- 
nary would lead to much higher values of ihb (see Hatzes et al. 
2003, Haghighipour2006). Although both fits correspond to sig- 
nificantly different dynamical systems, they yield practically the 
same residuals: Ot^)'''^ = 1.98. 

From the detection of exoplanets from RV data, we have 
learned (the hard way) that the best fit obtained with a limited 
data set does not necessarily correspond to the real configura- 
tion of the physical system. This is especially true for a small 
number of data points, or when they only cover a fraction of 
the orbital period. In other words, even if we specify a value for 
iriA, the real mass and orbit of niB does not necessarily have to 
coincide with the best fit. 

A way to estimate the possible range of solutions compatible 
with the observational data is to analyze the shape of the residual 
function for a series of orbital fits around the minimum value 
/fvmin- ^^^^ function shows a steep increase for small changes in 
the fitted parameters, then we may have a certain confidence that 
the best fit is probably very close to the real system. However, 
a shallow minimum could lead to a wide diversity of possible 
solutions with almost the same value of;^'^ and, consequently, to 
different configurations that are statistically indistinct. 

To test this idea for both choices of m^, we calculated a series 
of orbital solutions with predefined values of as between 12 and 
25 AU. Except for the pair {111^ ,aB), all the remaining parameters 
were free and determined with the minimization process. Results 
are shown in Figure |5]for niA = 1.59M0 and for iua = 1.1 8M0. 
The value of (xl)^^^ for each fit is shown in the top lefthand 
plot. Both families of solutions show similar minima, although 
displaced in semimajor axis. 

To estimate the limits of the Icr confidence level, we em- 
ployed the same procedure as employed in Giuppone et al. 
(2009). Given a best-fit algorithm with v degrees of freedom. 
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Fig. 5. Several multikeplerian orbital fits for y-Cephei consider- 
ing different values of niA and ub- Black lines show results for 
niA - I.59M0, while red curves assume itia = 1.1 8M0. The 
dashed horizontal line in the upper-left panel corresponds to the 
Icr confidence level around the global minimum of;^'^. 

the value of (xl)^^^ associated to the Icr confidence level is ap- 
proximately given by 



where (vlV^^ 



1 + A — 



(21) 



is the minimum value. Since our problem contains 
230 data points and 10 free parameters, we have v = 220. For 
- 2.0, equation (EB then gives (xl)\'^ - 2.1. 
Although both values only differ in ~ 5%, the range of 
possible solutions with (xl)^^^ < (xl)uj is surprisingly wide. 
Assuming the lower value for niA, the semimajor axis may be 
as small as 15 AU or as large as 22 AU. For niB - 1.59M0, the 
range is ub e [17, 23] AU. 

The top righthand plot of Figure |5]presents the values of 
that give the best fit for each value of gb. As expected, when 
the components are more separated, the mass of the binary in- 
creases to produce same magnitude in radial velocity. The lower 
lefthand plot gives the binary eccentricity e^. Larger semimajor 
axes are accompanied by more elliptic orbits. The curves show a 
rough resemblance to the locus of constant pericentric distance 
Pb, here drawn for one particular orbital solution. Since most of 
the RV data points are located near the pericenter of the binary's 
orbit (see Figure |4), the value of pb is much better defined than 
either as or es. 

Finally, the lower-right panel shows that the offset from the 
CfA data varies significantly as a function of the semimajor axis 
of binary. This behavior is not noted for the other offsets, which 
remain almost constant. All the best fit solutions showed almost 
no variation in the longitude of pericenter nor in the time of pas- 
sage through the periastron. 

Table|2]summarizes the minimum/maximum possible values 
of niB, Ob, and leading to orbital fits with residuals within 
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Table 2. Icr Confidence limits for y-Cephei-B 





triA = I.I8M0 


niA = 1.59Mo 




[0.25, 0.34] 
[15,22] 
[0.32, 0.48] 


[0.3,0.4] 
[17,23] 
[0.33,0.46] 



the Icr confidence level. It is important to stress that statistically 
all are equally compatible with the observational data. In the next 
section we attempt to elucidate whether these uncertainties are 
important to the dynamics of small planetesimals orbiting the 
primary star and, consequently, whether they could have an ef- 
fect on the planetary formation process. 

4. Secular dynamics of planetesimals in y-Cephei 

4. 1 . The forced eccentricity 

Planetary accretion requires low relative velocities which, in 
turn, implies similar orbits between colliding bodies. This con- 
dition may be satisfied if the orbital eccentricities are: (i) very 
small or (ii) very similar and the orbits are aligned. For small 
planetesimals where mutual perturbations are not crucial, the or- 
bital eccentricities are determined by a complex interplay be- 
tween several phenomena, including gas drag, collisional his- 
tory, and the gravitational effects of the secondary star (e.g. 
IVIarzari and Scholl 2000, Thebault et al. 2006, 2008). In the sec- 
ular approximation, these effects appear through the magnitude 
of the forced eccentricity ej. 

Thus, one way to study planetary accretion under different 
orbital solutions for y-Cephei would be to analyze the sensitivity 
of Cf to the set {mA,mB, 03,^3) compatible with the RV data. 
Results are presented in Figure |6] for two values of niA- Each 
panel shows level curves of constant forced eccentricity, as a 
function of the semimajor axis a of the planetesimal (abscissa), 
and for different semimajor axes a b for the binary pair (ordinate). 

Contrary to our expectations, the range of eccentricities ap- 
pears insensitive to the binary configuration. In all cases the val- 
ues of Cf extend from ~ 0.03 for the small semimajor axis to 
~ 0.077 for fl ^ 4 AU. Adopting a lower value of niA seems 
to cause a slight reduction in the interval, but the change is not 
very significant. Consequently, and at least from this preliminary 
analysis, there appears to be no indication that different config- 
urations for y-Cephei could cause major changes in the secular 
dynamics of small bodies and, therefore, on the accretional pos- 
sibilities of a planetesimal swarm. 

With hindsight, perhaps this result is not at all unexpected. 
Since all orbital solutions for y-Cephei lead to practically the 
same amplitude in the RV signal, it is understandable that dif- 
ferent values for the set (m^i, m^, a^, e^) should also generate 
similar perturbative effects on other hypothetical bodies in the 
system; e.g. small planetesimals orbiting m^. 

4.2. Simulations of individuai particles with gas drag 

A better test for the effects of different orbital fits on the secu- 
lar dynamics is to compare the evolution of small planetesimals 
under the effects of a nonlinear drag force from a circumstellar 
gas disk centered on niA- We employ the same expression for 
the dissipative force as discussed in Beauge et al. (2010). For 
the gaseous disk we assume a linear surface density profile with 
a total mass of 3 Jupiter masses, and an outer edge located at 




1 1.5 2 2.5 3 3.5 4 

a [AU] 



h 


'1 

1.18 \ 


\ 




' 1 


\ 


- \ 0.04 


■V-v- 

\ 0.05 

\v 




0.07 
,\ 1 


0.073 ] ^ 

, r — r 




1.5 


2 


2.5 


3 


3.5 ^ 



a [AU] 
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Figure|5] Values of as for each best fit are shown with horizontal 
dashed lines. 



5 AU. We consider planetesimals with a volumetric density of 
p = 3 g/cml 

Figure |7] shows the result of numerical simulations of four 
different test planetesimals, with radius between s - 0.5 km 
(top) and s - 10 km (bottom). All were initially placed in cir- 
cular orbits. The initial semimajor axis was equal to a = 4 AU 
for the first three integrations, and a = 3 AU for the largest plan- 
etesimal. We assumed a gas disk with constant eccentricity of 
gg = 0.05 and a rigid-body retrograde precession with a period 
onOOOyrs. 

We considered two different orbital fits for the stellar com- 
ponents. The mass and orbital parameters for the secondary star 
were taken from equations (l20t which are the best fit solutions 
for each value of m^. 

The lefthand panels show the orbital eccentricity as a func- 
tion of the semimajor axis. The orange curves indicate the forced 
eccentricity e/, as obtained from our second-order model, for 
each value of niA ■ The righthand plots show the evolution of the 
longitude of pericenter m. The origin was shifted so that -utb = 0. 

We note the existence of three distinct regions in the semi- 
major axis domain. For a <2 AU, the planetesimals show short- 
period oscillations around equilibrium values ggq and tzTeq. These 
equilibrium values depend on size. Small planetesimals show 
^eq < and tiTeq < 0. Larger bodies, however, appear coupled 
with the conservative secular solution (k, h) - {ej, 0). 
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Fig. 7. Orbital evolution of four different size planetesimals un- 
der the effects of a nonlinear gas drag in the y-Cephei system. 
Black dots correspond to iua = 1.59Mq and gray to niA - 
1.1 SMq. In both cases the parameters of mg are those given by 
the best fits and detailed in equations ( |20] |. The gaseous disk has 
a constant eccentricity of e/ = 0.05 and a rigid retrograde pre- 
cession rate with period 2;7r/|^g| = 1000 yr In the lefthand pan- 
els, the orange curve shows the forced eccentricities as a function 
of the semimajor axis. In the righthand plots, the orange curves 
marks m - Q. 



The second region lies roughly between 2 AU and 3 AU, and 
is characterized by very similar equilibrium values of both e and 
m for all planetesimals with radius i > 0. 1 km. Moreover, the 
equilibrium eccentricities are virtually indistinguishable from 
the forced eccentricity e /. If planetary accretion is possible with 
these disk parameters, this region appears to be the most promis- 
ing since orbital dispersion is kept at a minimum. The limit be- 
tween both regions (here at a ~ 2 AU) depends on the surface 
density profile adopted for the disk. In our simulations we used a 
linear dependence with the radial distance (Beauge et al. 2010), 
which translates to high densities close to the central star and low 
values beyond 3 AU (see Figure 3 of Kley and Nelson 2008). 

A third region is located beyond ~ 3 AU. Although the sec- 
ular dynamics also appear similar for all values of s, the sim- 
ulations show large-amplitude oscillations. These are not only 
caused by short-period terms but also from several high-order 
mean-motion resonances between the particles and mg. In the 
plots these commensurabilities can be seen as spikes where the 
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Fig. 8. Example of trapping of a s = 10 km planetesimal in 
a 10/1 mean-motion resonance with the secondary star of y- 
Cephei, due to a nonlinear gas drag. Although the orbital migra- 
tion is divergent, capture still occurs and leads to an apparently 
stable configuration. The resonant angle is = IQAb - A- 9mB. 



eccentricity is temporarily excited. Without a detailed analysis, 
it is not possible to establish whether these commensurabilities 
will inhibit or favor accretion. Although most of our simulations 
have shown scattering effects and significant orbital misalign- 
ment between resonant and non-resonant orbits, we have also 
found cases of resonant trapping. This appears to be a high- 
probability outcome for s >5 km. 

An example is shown in Figure[8]for a 10 km body placed in 
an initial circular orbit with a - A AU. After an initial decay in 
the semimajor axis, the body is captured in a 10/1 mean-motion 
resonance with the binary component. Both the resonant angle 



e = \QAb -a- 9mB 



(22) 



and the difference in longitudes of pericenter librate around zero, 
although there seems to be a slow departure towards asymmet- 
ric librations at the end of the simulation. The resonant solution 
seems very stable, at least for timescales between 10^ and 10^ 
years. We have found similar outcomes in other commensura- 
bilities, such as the 11/1 and 12/1, always for slow dissipative 
effects (i.e. large bodies). What is curious about this result is that 
trapping occurs during a divergent migration; in other words, 
the non-conservative force increases the separation of the bodies 
involved. Classical works (e.g. Neishtadt 1975, Henrard 1982, 
Beauge and Ferraz-lVIello 1993, Nelson and Papaloizou 2002) 
predict that capture is only possible in cases of convergent mi- 
gration and, thus, the behavior shown in Figure |8] should not oc- 
cur 

The only reference we have been able to find describing sim- 
ilar findings is an abstract of a presentation in the 2007 DDA 
meeting (Hamilton and Zhang 2007). Although no details are 
available, it appears that divergent trapping is possible in high- 
order mean-motion resonances and with high-eccentricity per- 
turbers. This may explain why the librating resonant angle in- 
cludes the longitude of pericenter of the perturber instead of the 
planetesimal. However, a deeper analysis is necessary before we 
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Fig. 9. Same as Figure |7] but assuming a static gas disk anti- 
aligned with the secondary star niB- 

are able to understand this phenomena and establish its possible 
importance in planetary formation. 

The results shown in Figure Qfor both binary configurations 
show almost identical results. Although the forced eccentricity 
is slightly higher for niA - 1.1 8M0 for most of the semimajor 
axis domain, the same three regions exist in both cases, and it is 
plausible to assume that the collisional evolution of a planetesi- 
mal swarm should be similar Figure |9]shows a second series of 
integrations with the same four planetesimals as before, but this 
time we considered no precession for the gaseous disk. We also 
assumed that TiJg - vjb - n, i.e. the disk is anti-aligned with the 
binary companion. This is consistent with the hydrosimulations 
of Marzari et al. (2009b) for disks with significant self-gravity. 
Although the magnitude of the resonant and short-period oscil- 
lations appear larger for ma - 1.1 8M0, the averaged secular 
dynamics is similar in both cases, and the same three regions 
discussed previously are also present. 

5. Conclusions 

In this paper we have presented a second-order theory for the 
secular dynamics of massless particles orbiting a central star and 
perturbed by a secondary stellar component with high eccentric- 
ity. Only coplanar motion is considered. This dynamical prob- 
lem is applicable to the motion of small planetesimals in tight 
binary systems such as y-Cephei. Although the resulting expres- 



sions for the forced eccentricity e f and secular frequency g are 
complex, we were able to extend the empirical approximation 
originally introduced by Thebault et al. (2006) and deduce sim- 
ple analytical formulas for both quantities. 

The forced eccentricity, in particular, shows significant dif- 
ferences to the classical model (e.g. Heppenheimer 1978). While 
the first-order equations predict a linear dependence with the 
semimajor axis, numerical simulations show a quadratic func- 
tional form, and e j may actually reach a plateau for high values 
of a. Our model reproduces this behavior with good precision. 

We also analyzed the reliability of the best fits presented in 
the literature for the stellar components of y-Cephei. We found 
that the best solution depends on the adopted mass for the cen- 
tral star, and even for a fixed value of niA there may be many dif- 
ferent configurations compatible with the observations. This is 
expected, since the radial velocity data covers less than one or- 
bital period of the system. However, we have also found that the 
dynamics of small planetesimals appears to be only weakly de- 
pendent on the particular solution adopted for the binary. A com- 
parative study of the evolution of planetesimals under the effects 
of gas drag from a circumstellar {niA centered) gas disk shows 
similar evolutionary paths. Although our integrations have not 
covered many initial conditions or disk parameters, we suspect 
that the accretional process of a planetesimal swarm should be 
practically equivalent in any case. Consequently, we believe that 
the difficulties in explaining planetary formation in tight binary 
systems cannot be attributed to uncertainties in the orbital fits. 
The solution must be found elsewhere, and the search is pre- 
cisely what makes this problem intriguing. 

Finally, we also presented a curious case of resonant trapping 
in divergent migration. There is practically no reference to this 
behavior in the literature and, although its importance in plane- 
tary formation is difficult to evaluate, we nevertheless believe the 
phenomena is intrinsically interesting and merits further analy- 
sis. 
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